# -*- coding: utf-8 -*-

"""
SAL方法
@author: qzhang
"""

import pymongo
import numpy as np
import string
import datetime
import flask
from profilehooks import profile

# 计算场相似
# 传入预报结果站号、经纬度信息、预报值
# 预报结果以数组形式存储
# fore=[站号,经度,纬度,预报值]
# 传入观测结果站号、经纬度信息、观测值
# 观测结果以数组形式存储
# obs=[站号,经度,纬度,预报值]
# 传入需要计算的变量的量级
# 量级为int型
# 例如降水的小雨、中雨、大雨
@profile
def field_corr(fore,obs):
    corrcoef_matrix=np.corrcoef(fore[:,3],obs[:,3])#用numpy现有模块计算自相关与交叉相关系数
	field_error=corrcoef_matrix[0,1]#corrcoef_matrix[0,1]就是场相关系数
	return field_error

# 计算面积重合比(S)
# 传入预报结果站号、经纬度信息、预报值
# 预报结果以数组形式存储
# fore=[站号,经度,纬度,预报值]
# 传入观测结果站号、经纬度信息、观测值
# 观测结果以数组形式存储
# obs=[站号,经度,纬度,预报值]
# 传入需要计算的变量的量级
# 量级为int型
# 例如降水的小雨、中雨、大雨
@profile
def area_overlap(fore,obs,level):
	if level==0:#小雨级别雨量等级
		thresh=[0,10]
	if level==1:#中雨级别雨量等级
		thresh=[10,25]
	if level==2:#大雨级别雨量等级
		thresh=[25,50]
	if level==3:#暴雨级别雨量等级
		thresh=[50,999]
	count_fore=0#预报值在某一级别内的站点个数
	count_obs=0#实况值在某一级别内的站点个数
	count_right=0#报准的站点个数
	count=0
	for cell in fore:#按照从数据库取出来的结果遍历
		flag_obs=False#如果站点雨量在某一量级之内为True，反之则为False
		flag_fore=False#如果站点雨量在某一量级之内为True，反之则为False
		if cell[3]<thresh[1] and cell[3]>=thresh[0]:#判断某一站点预报值是否在范围内
			count_fore+=1
			flag_fore=True
		else:
			flag_fore=False
		if obs[count,3]<thresh[1] and obs[count,3]>=thresh[0]:#判断某一站点实况值是否在范围内
			count_obs+=1
			flag_obs=True
		else:
			flag_obs=False
		if flag_obs==True and flag_fore==True:#判断站点是否报准
			count_right+=1
		count+=1
	if count_obs!=0 and count_fore!=0:
		overlap=float(count_right)/float(count_obs+count_fore)#计算准确率
	else:
		overlap=None
	return overlap
	
# 计算强度偏差(A)
# 传入预报结果站号、经纬度信息、预报值
# 预报结果以数组形式存储
# fore=[站号,经度,纬度,预报值]
# 传入观测结果站号、经纬度信息、观测值
# 观测结果以数组形式存储
# obs=[站号,经度,纬度,预报值]
# 传入需要计算的变量的量级
# 量级为int型
# 例如降水的小雨、中雨、大雨
@profile
def strength_error(fore,obs,level):
	if level==0:#小雨级别雨量等级
		thresh=[0,10]
	if level==1:#中雨级别雨量等级
		thresh=[10,25]
	if level==2:#大雨级别雨量等级
		thresh=[25,50]
	if level==3:#暴雨级别雨量等级
		thresh=[50,999]
	count_fore=0#预报结果在某一量级之内的站点个数
	count_obs=0#实况值在某一量级内的站点个数
	sum_fore=0#预报值在某量级内的预报值总和
	sum_obs=0#观测值在某量级内的预报值总和
	count=0
	for cell in fore:#按照从数据库取出来的结果遍历
		if cell[3]<thresh[1] and cell[3]>=thresh[0]:#判断某一站点预报值是否在范围内
			count_fore+=1
			sum_fore+=cell[3]
		if obs[count,3]<thresh[1] and obs[count,3]>=thresh[0]:#判断某一站点实况值是否在范围内
			count_obs+=1
			sum_obs+=obs[count,3]
		count+=1
	if count_fore==0 and count_obs!=0:
		ave_stren_err=0-(sum_obs/count_obs)
	if count_fore!=0 and count_obs==0:
		ave_stren_err=(sum_fore/count_fore)-0
	if count_fore!=0 and count_obs!=0:
		ave_stren_err=(sum_fore/count_fore)-(sum_obs/count_obs)#平均强度误差=预报平均值-实况平均值
	if count_fore==0 and count_obs==0:
		ave_stren_err=None
	max_stren_err=np.max(fore[:,3])-np.max(obs[:,3])#极大值强度误差=预报极大值-实况极大值
	min_stren_err=np.min(fore[:,3])-np.min(obs[:,3])#极小值强度误差=预报极小值-实况极小值
	return [ave_stren_err,max_stren_err,min_stren_err]
		
# 计算质心偏差(A)
# 传入预报结果站号、经纬度信息、预报值
# 预报结果以数组形式存储
# fore=[站号,经度,纬度,预报值]
# 传入观测结果站号、经纬度信息、观测值
# 观测结果以数组形式存储
# obs=[站号,经度,纬度,预报值]
# 传入需要计算的变量的量级
# 量级为int型
# 例如降水的小雨、中雨、大雨
@profile
def location_error(fore,obs,level):
	if level==0:#小雨级别雨量等级
		thresh=[0,10]
	if level==1:#中雨级别雨量等级
		thresh=[10,25]
	if level==2:#大雨级别雨量等级
		thresh=[25,50]
	if level==3:#暴雨级别雨量等级
		thresh=[50,999]
	mass_fore=0#预报值在某一量级之内的总和
	mass_obs=0#实况值在某一量级之内的总和
	lat_fore=0#预报值在某一量级之内的所有站点纬度降水权重总和
	lat_obs=0#实况值在某一量级之内的所有站点纬度降水权重总和
	lon_fore=0#预报值在某一量级之内的所有站点纬度降水权重总和
	lon_obs=0#实况值在某一量级之内的所有站点纬度降水权重总和
	count=0
	for cell in fore:#按照从数据库取出来的结果遍历
		if cell[3]<thresh[1] and cell[3]>=thresh[0]:#判断某一站点预报值是否在范围内
			mass_fore+=cell[3]
			lon_fore+=cell[1]*cell[3]
			lat_fore+=cell[2]*cell[3]
		if obs[count,3]<thresh[1] and obs[count,3]>=thresh[0]:#判断某一站点实况值是否在范围内
			mass_obs+=obs[count,3]
			lon_obs+=obs[count,1]*obs[count,3]
			lat_obs+=obs[count,2]*obs[count,3]
		count+=1
	if mass_fore!=0 and mass_obs!=0:
		lat_error=(lat_fore/mass_fore)-(lat_obs/mass_obs)#质心纬度偏差=预报质心纬度-实况质心纬度
		lon_error=(lon_fore/mass_fore)-(lon_obs/mass_obs)#质心经度偏差=预报质心经度-实况质心经度
	if mass_fore==0 and mass_obs==0:
		lat_error=None
		lon_error=None
	if mass_fore!=0 and mass_obs==0:
		lat_error=(lat_fore/mass_fore)-0
		lon_error=(lon_fore/mass_fore)-0
	if mass_fore==0 and mass_obs!=0:
		lat_error=0-(lat_obs/mass_obs)
		lon_error=0-(lon_obs/mass_obs)
	return [lon_error,lat_error]

#不同类型的站点列表配置文件路径
station_path="/home/qzhang/桌面/station.ini"

#服务器配置文件
ini_path="/home/qzhang/桌面/dunoinfo.ini"

#前端传回的请求
get={"fdate":"20150101",								#		string	开始时间
			"send":"H16",											#		list			发报时间
			"timestep":"00_24",								#		string	预报时间步长
			"ftype":"FD",											#		string	预报类型
			"var":"Prcp"}											#		string	查询变量名称缩写		("Prcp":降水,"WDir":风向,"Wspd":风速)

#连接服务器并获取数据
@profile
def get_value(get,ini_path,station_path):
	"""读取站点配置文件"""
	flag_station=open(station_path,'r')
	station_data=flag_station.read().split("\n")
	station_num=len(station_data)
	data_fore=np.zeros([station_num,4])
	data_obs=np.zeros([station_num,4])
	"""读配置文件并连接服务器"""
	flag_ini=open(ini_path,"r")
	ini_data=flag_ini.read().split("\n")
	usrnm=ini_data[0]
	pswd=ini_data[1]
	dm=ini_data[2]
	dbnm=ini_data[3]
	classnm=ini_data[4]
	print "database info :",usrnm,pswd,dm,dbnm,classnm
	"""连接服务器"""
	connection=pymongo.MongoClient(dm)
	if connection[dbnm].authenticate(usrnm,pswd)==True:
		db=connection[dbnm]
	else:
		print "Connect To Sever Error, Check  ini file"
	"""从服务器读数据"""
	count=0
	for cell in station_data:
		station_info=cell.split("\t")
		if station_info==['']:
			break
		else:
			data_fore[count,0]=int(string.atof(station_info[0][1:]))
			data_fore[count,1]=string.atof(station_info[1])/10000
			data_fore[count,2]=string.atof(station_info[2])/10000
			data_obs[count,:]=data_fore[count,:]
			data=db[classnm].find_one({"AWSID":station_info[0],"FDate":get["fdate"],"Send":get["send"],"FVali":get["ftype"]+get["timestep"]})
			if data==None:
				count+=1
				continue
			else:
				dict_var="Sm_"+get["var"]
				dict_fore="Fore_"+get["var"]
				dict_obs="Ob_"+get["var"]
				data_fore[count,3]=data[dict_var][dict_fore]
				data_obs[count,3]=data[dict_var][dict_obs]
				count+=1
	return data_fore,data_obs
	
#test
level=1
fore,obs=get_value(get,ini_path,station_path)
print field_corr(fore,obs)
print area_overlap(fore,obs,level)
print strength_error(fore,obs,level)
print location_error(fore,obs,level)

		
